program hmap

 use mod_utils, only: &
   t_realtime, my_second

 use mod_messages, only: &
   mod_messages_constructor, &
   mod_messages_destructor,  &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_constructor, &
   mod_kinds_destructor,  &
   wp

 use mod_fu_manager, only: &
   mod_fu_manager_constructor, &
   mod_fu_manager_destructor,  &
   new_file_unit

 use mod_octave_io, only: &
   mod_octave_io_constructor, &
   mod_octave_io_destructor,  &
   real_format,    &
   write_octave,   &
   read_octave,    &
   read_octave_al, &
   locate_var
 
 use mod_linal, only: &
   mod_linal_constructor, &
   mod_linal_destructor

 use mod_perms, only: &
   mod_perms_constructor, &
   mod_perms_destructor

 use mod_octave_io_perms, only: &
   mod_octave_io_perms_constructor, &
   mod_octave_io_perms_destructor

 use mod_numquad, only: &
   mod_numquad_constructor, &
   mod_numquad_destructor,  &
   get1d_gauss_nodes, get1d_gausslobatto_nodes

 use mod_sparse, only: &
   mod_sparse_constructor, &
   mod_sparse_destructor

 use mod_octave_io_sparse, only: &
   mod_octave_io_sparse_constructor, &
   mod_octave_io_sparse_destructor

  use mod_sympoly, only: &
   mod_sympoly_constructor, &
   mod_sympoly_destructor

 use mod_octave_io_sympoly, only: &
   mod_octave_io_sympoly_constructor, &
   mod_octave_io_sympoly_destructor

 use mod_mpi_utils, only: &
   mod_mpi_utils_constructor, &
   mod_mpi_utils_destructor,  &
   mpi_comm_world, wp_mpi, &
   mpi_init, mpi_finalize, &
   mpi_comm_size, mpi_comm_rank

 use mod_master_el, only: &
   mod_master_el_constructor, &
   mod_master_el_destructor,  &
   mod_master_el_initialized

 use mod_grid, only: &
   mod_grid_constructor, &
   mod_grid_destructor,  &
   mod_grid_initialized, &
   t_grid, t_ddc_grid,   &
   new_grid, clear,      &
   affmap,               &
   write_octave
 
 use mod_bcs, only: &
   mod_bcs_constructor, &
   mod_bcs_destructor,  &
   mod_bcs_initialized, &
   t_bcs,                     &
   t_b_v,   t_b_s,   t_b_e,   &
   new_bcs, clear,            &
   write_octave

 use mod_base, only: &
   mod_base_constructor, &
   mod_base_destructor,  &
   mod_base_initialized, &
   t_base, clear, &
   write_octave
 
 use mod_state_vars, only: &
   mod_state_vars_constructor, &
   mod_state_vars_destructor,  &
   mod_state_vars_initialized, &
   c_stv

 use mod_hmap_state, only: &
   mod_hmap_state_constructor, &
   mod_hmap_state_destructor,  &
   mod_hmap_state_initialized, &
   t_hmap_state, clear, &
   new_hmap_state,      &
   write_octave
 
 use mod_output_control, only: &
   mod_output_control_constructor, &
   mod_output_control_destructor

 use mod_time_integrators, only: &
   mod_time_integrators_constructor, &
   mod_time_integrators_destructor,  &
   mod_time_integrators_initialized, &
   c_ode, t_ti_init_data, t_ti_step_diag,     &
   i_ti_init,    i_ti_step,    i_ti_clean,    &
   thetam_init,  thetam_step,  thetam_clean,  &
   bdf1_init,    bdf1_step,    bdf1_clean,    &
   bdf2_init,    bdf2_step,    bdf2_clean,    &
   bdf3_init,    bdf3_step,    bdf3_clean,    &
   erb2_init,    erb2_step,    erb2_clean
 
 use mod_hmap_ode, only: &
   mod_hmap_ode_constructor, &
   mod_hmap_ode_destructor,  &
   mod_hmap_ode_initialized, &
   t_hmap_ode
   
 use mod_linsolver, only: &
   mod_linsolver_constructor, &
   mod_linsolver_destructor

!-----------------------------------------------------------------------

 implicit none

!-----------------------------------------------------------------------

 ! Parameters
 integer, parameter :: max_char_len = 1000
 character(len=*), parameter :: input_file_name = 'hmap.in'
 character(len=*), parameter :: this_prog_name = 'hmap'

 ! Main variables
 integer :: mpi_id, mpi_nd
 logical :: write_grid, write_sys, compute_error_norms
 integer :: k, nbound
 character(len=max_char_len) :: grid_file, in_testname, in_basename
 type(t_base) :: base
 type(t_grid), target :: grid
 type(t_ddc_grid) :: ddc_grid
 type(t_bcs), target :: bcs

 ! Time stepping
 logical :: l_checkpoint, l_output, l_statstics
 integer :: nstep, n, n_out
 real(wp) :: tt_sta, tt_end, dt, t_n, t_nm1, time_last_out, dt_out, &
   time_last_stats, dt_stats, time_last_check, dt_check

 ! ODE variables and state vectors
 type(t_hmap_state) :: uuu0, uuun
 type(t_hmap_ode)  :: hmap_ode
 type(t_ti_init_data) :: ti_init_data
 type(t_ti_step_diag) :: ti_step_diag
 procedure(i_ti_init),  pointer :: ti_init  => null()
 procedure(i_ti_step),  pointer :: ti_step  => null()
 procedure(i_ti_clean), pointer :: ti_clean => null()
 logical :: l_output1, l_output2

 ! Auxiliary variables
 integer :: fu, ierr
 real(t_realtime) :: t00, t0, t1
 character(len=10*max_char_len) :: message(3)


  call mod_messages_constructor()

  call mod_kinds_constructor()
  
  call mod_fu_manager_constructor()

  call mod_octave_io_constructor()

  call mod_linal_constructor()

  call mod_perms_constructor()
  call mod_octave_io_perms_constructor()

  call mod_numquad_constructor()

  call mod_sparse_constructor()
  call mod_octave_io_sparse_constructor()

  call mod_sympoly_constructor()
  call mod_octave_io_sympoly_constructor()

  call mpi_init(ierr)
  call mod_mpi_utils_constructor()
  call mpi_comm_size(mpi_comm_world,mpi_nd,ierr)
  call mpi_comm_rank(mpi_comm_world,mpi_id,ierr)

  call mod_master_el_constructor()

  call mod_grid_constructor()

  call mod_bcs_constructor()

  call mod_base_constructor()

  call mod_state_vars_constructor()

  call mod_hmap_state_constructor()

  call mod_output_control_constructor('basename')

  call mod_time_integrators_constructor()

  call mod_linsolver_constructor()

!  !---------------------------------------------------------------------
!  ! Define the grid
!
!  ! Setup for the Landau damping
!  ! 1D
!  !allocate( xx(1) , vv(1) )
!  !allocate( xx(1)%x(51) , vv(1)%x(81) )
!  !xx(1)%x = (/ ( 4.0_wp*pi/50.0_wp*real(i,wp), i=0,50) /)
!  !vv(1)%x = (/ ( 20.0_wp/80.0_wp*real(i,wp)-10.0_wp, i=0,80) /)
!  ! 2D
!  allocate( xx(2) , vv(2) , nex(2) , nev(2) )
!  nex = (/8,4/) ; nev = (/12,8/)
!  allocate( xx(1)%x(nex(1)+1) , xx(2)%x(nex(2)+1) , &
!            vv(1)%x(nev(1)+1) , vv(2)%x(nev(2)+1) )
!  xx(1)%x = (/ ( 4.0_wp*pi/real(nex(1),wp)*real(i,wp), i=0,nex(1)) /)
!  xx(2)%x = (/ ( 4.0_wp*pi/real(nex(2),wp)*real(i,wp), i=0,nex(2)) /)
!  !vv(1)%x = (/( 18.0_wp/real(nev(1),wp)*real(i,wp)-9.0_wp, i=0,nev(1))/)
!  !vv(2)%x = (/( 18.0_wp/real(nev(2),wp)*real(i,wp)-9.0_wp, i=0,nev(2))/)
!  ! Adapted v grid
!  !vv(1)%x = (/ -1.0_wp , -0.5_wp , -0.25_wp , -0.125_wp , 0.0_wp , &
!  !              0.125_wp ,  0.25_wp ,  0.5_wp ,  1.0_wp /) * 9.0_wp
!  !vv(2)%x = vv(1)%x
!  ! Optimized grid
!  vv(1)%x = 9.0_wp * (/                                              &
!    -1.0_wp , -0.5_wp , -0.3_wp , -0.225_wp , -0.15_wp , -0.075_wp , &
!  0.0_wp,0.075_wp,0.15_wp ,  0.225_wp ,  0.3_wp  , 0.5_wp ,  1.0_wp /)
!  vv(2)%x = (/ -1.0_wp , -0.5_wp , -0.25_wp , -0.125_wp , 0.0_wp , &
!                0.125_wp ,  0.25_wp ,  0.5_wp ,  1.0_wp /) * 9.0_wp
!
!  ! Setup for the "bump on tail" testcase
!  !allocate( xx(1) , vv(1) , nex(1) , nev(1) )
!  !nex(1) = 15; nev = 16
!  !allocate( xx(1)%x(nex(1)+1) , vv(1)%x(nev(1)+1) )
!  !xx(1)%x = (/ ( 20.0_wp*pi/real(nex(1),wp)*real(i,wp), i=0,nex(1)) /)
!  !vv(1)%x = (/(18.0_wp/real(nev(1),wp)*real(i,wp)-9.0_wp, i=0,nev(1))/)
!
!  call new_grid(grid,xx,vv)
!
!  call new_file_unit(fu,ierr)
!  open(fu,file='test-setup.octxt',                                     &
!       status='replace',action='write',form='formatted',iostat=ierr)
!  call write_octave(grid,'grid',fu)
!  !---------------------------------------------------------------------
!
!  !---------------------------------------------------------------------
!  ! Define the basis
!  alpha_rec = .true.
!  call new_base(base(1),5,grid%d,6)
!  call new_base(base(2),7,grid%d)
!
!  call fill_xv_points()
!  ! Use these fields in octave as
!  !   [X,Y]=meshgrid(x_points(:),v_points(:));
!  !   Z=uuu.f(:)(xv_idx);
!  !   surf(X,Y,Z')
!  !
!  ! And to make plots in 2D:
!  ! plot3(x_points(1,:,:)(:),x_points(2,:,:)(:),uuu.r(:),'o')
!  ! quiver(x_points(1,:,:)(:),x_points(2,:,:)(:), ...
!  !        uuu.et(1,1,:,:)(:),uuu.et(2,1,:,:)(:), 0)
!  ! To plot the distribution function in v
!  ! figure
!  ! hold on
!  ! iex = 1; ix = 1
!  ! pk = v_base.pk
!  ! for ie=1:grid.ne
!  !   if(grid.e(ie).e1(1)==iex)
!  !     f = reshape(uuu.f(ix,:,ie),[pk pk])
!  !     x = reshape(v_points(1,:,grid.e(ie).e1(2)),[pk pk])
!  !     y = reshape(v_points(2,:,grid.e(ie).e1(2)),[pk pk])
!  !     surf(x,y,f)
!  !   end
!  ! end
!  call write_octave(alpha_rec,'alpha_rec',fu)
!  call write_octave(base(1)  ,'x_base'   ,fu)
!  call write_octave(base(2)  ,'v_base'   ,fu)
!  call write_octave(xv_idx   ,'xv_idx'   ,fu)
!  call write_octave(x_points ,'x_points' ,fu)
!  call write_octave(v_points ,'v_points' ,fu)
!  call write_octave(xv_points,'xv_points',fu)
!  call write_octave(xx_points,'xx_points',fu)
!  !---------------------------------------------------------------------
!
!  !---------------------------------------------------------------------
!  ! Linear solver
!  !---------------------------------------------------------------------
!
  !---------------------------------------------------------------------
  ! Allocations and initial condition
  call mod_hmap_ode_constructor()
 
  call new_hmap_state(uuu0,grid,base)
  call new_hmap_state(uuun,grid,base)
 
  uuu0%d(:,1) = (/1.0_wp,0.0_wp/); uuu0%q = 0.0_wp
  !---------------------------------------------------------------------

  !---------------------------------------------------------------------
  ! Time integration
  !
  !   0        1                 n             nstep
  !   |--------|--------|--------|--------|------|
  ! tt_sta            t_nm1     t_n            tt_end
  !                     | step_n |
  !
  tt_sta = 0.0_wp; tt_end = 100.0_wp; dt = 1.0_wp
  nstep = ceiling(tt_end/dt)
  n_out = 0
  call write_out1(0,tt_sta,n_out,uuu0,ti_step_diag)
  call write_out2(tt_sta,uuu0,.true.)
  call write_out2(tt_sta,uuu0,.false.)

  ti_init  => bdf3_init
  ti_step  => bdf3_step
  ti_clean => bdf3_clean
  ti_init_data%mpi_id = mpi_id
  ti_init_data%mpi_comm = mpi_comm_world
  call ti_init(hmap_ode,dt,tt_sta,uuu0,ti_init_data,ti_step_diag)
  time_loop: do n=1,nstep

    t_nm1 = tt_sta + real(n-1,wp)*dt
    t_n   = t_nm1 + dt

    call ti_step(hmap_ode,uuun,t_nm1,uuu0,ti_step_diag)
    call uuu0%copy(uuun)

    l_output1 = (mod(n,1).eq.0).or.(n.eq.1).or.(n.eq.nstep)
    if(l_output1) then
      write(*,*) 'Step ',n,' of ',nstep
      n_out = n_out + 1
      call write_out1(n,t_n,n_out,uuu0,ti_step_diag)
    endif

    l_output2 = (mod(n,1).eq.0).or.(n.eq.1).or.(n.eq.nstep)
    if(l_output2) call write_out2(t_n,uuu0,.false.)

  enddo time_loop
  call ti_clean(hmap_ode,ti_step_diag)
  call clear(uuu0)
  call clear(uuun)
  !---------------------------------------------------------------------


  ! ?? call clear(grid)

  call mod_hmap_ode_destructor()

  call mod_linsolver_destructor()

  call mod_time_integrators_destructor()

  call mod_output_control_destructor()

  call mod_hmap_state_destructor()

  call mod_state_vars_destructor()

  call mod_base_destructor()

  call mod_bcs_destructor()

  call mod_grid_destructor()

  call mod_master_el_destructor()

  call mod_mpi_utils_destructor()
  call mpi_finalize(ierr)

  call mod_octave_io_sympoly_destructor()
  call mod_sympoly_destructor()

  call mod_octave_io_sparse_destructor()
  call mod_sparse_destructor()

  call mod_numquad_destructor()

  call mod_octave_io_perms_destructor()
  call mod_perms_destructor()

  call mod_linal_destructor()

  call mod_octave_io_destructor()

  call mod_fu_manager_destructor()

  call mod_kinds_destructor()

  call mod_messages_destructor()

contains


 subroutine write_out1(n,t,n_out,uuu,sd)
  integer, intent(in) :: n, n_out
  real(wp), intent(in) :: t
  type(t_hmap_state), intent(in) :: uuu
  type(t_ti_step_diag), intent(in) :: sd

  integer :: fu, ierr
  character(len=*), parameter :: time_stamp_format = '(i4.4)'
  character(len=4) :: time_stamp

   call new_file_unit(fu,ierr)
   write(time_stamp,time_stamp_format) n_out
   open(fu,file='output-'//time_stamp//'.octxt', &
        status='replace',action='write',form='formatted',iostat=ierr)
    call write_octave(n  ,'n'  ,fu)
    call write_octave(t  ,'t'  ,fu)
    call write_octave(uuu,'uuu',fu)
   close(fu,iostat=ierr)

 end subroutine write_out1

 !> Write "high-frequency" output
 subroutine write_out2(t,uuu,init)
  real(wp), intent(in) :: t
  type(t_hmap_state), intent(in) :: uuu
  logical, intent(in) :: init

  character(len=*), parameter :: outfile = 'fast-output.octxt'
  character(len=5) d_size

   if(init) then
     call new_file_unit(fu,ierr)
     open(fu,file=outfile, status='replace', &
       action='write',form='formatted')
     ! Nothing to write: we just make sure the file is empty
     close(fu)
     return
   endif

   write(d_size,'(i5)') 4
   call new_file_unit(fu,ierr)
   open(fu,file=outfile, status='old',                    &
     action='readwrite',form='formatted',position='append')
    write(fu,'('//d_size//'('//real_format//')'//')') (/t,uuu%d,uuu%q/)
   close(fu)

 end subroutine write_out2

end program hmap
